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EXECUTIVE  SUMMARY 


A  critical  step  in  modeling  groundwater  flow  and/or  contaminant  transport  at  a  site  is  to  obtain  an 
accurate  3-D  model  of  the  hydraulic  conductivity  of  the  subsurface.  Past  practice  has  often 
involved  assigning  average  values  of  hydraulic  conductivity  K  to  simple  layered  models  of  the 
subsurface.  However  it  now  widely  recognized  that  using  averaged  values  of  K  is  insufficient; 
the  spatial  variability  in  K,  the  natural  heterogeneity  of  the  subsurface,  has  a  significant  effect  on 
groundwater  flow  and  is  an  important  parameter  in  modeling  contaminant  transport  The  challenge 
is  to  develop  methods  that  can  be  used  to  accurately,  and  cost-effectively  characterize  the 
subsurface  in  terms  of  its  hydraulic  properties.  Information  is  needed  about  both  the  magnitude 
and  the  spatial  variability  in  properties.  We  have  investigated  the  use  of  measurement  of  dielectric 
properties  as  a  means  of  obtaining  both  these  types  of  information. 

A  central  goal  of  our  research  is  to  develop  an  understanding  of  the  link  between  dielectric 
properties  and  hydraulic  properties  of  the  subsurface.  Measurements  of  dielectric  properties  of  the 
subsurface  can  be  obtained  non-invasively  from  ground  penetrating  radar  data.  The  dielectric 
constant  of  a  volume  of  the  subsurface  controls  the  velocity  of  the  electromagnetic  wave;  changes 
in  the  dielectric  constant  of  the  subsurface  determine  the  magnitude  of  the  reflections  seen  in  GPR 
data.  In  this  way  a  "map"  of  the  dielectric  constant  of  the  subsurface  can  be  obtained  from  GPR 
data.  From  an  applied  perspective  our  research  can  be  described  as  addressing  the  question:  How 
can  we  use  these  measurements  of  dielectric  constant  to  characterize  hydraulic  conductivity  at  a 
site?  Our  approach  is  a  combination  of  field  experiments,  laboratory  experiments,  and  theoretical 
modeling. 

One  focus  of  our  research  has  been  to  work  with  GPR  data  and  to  develop  data  analysis  methods 
that  allow  us  to  extract  information  about  the  spatial  variability  in  hydraulic  properties.  Specifically 
we  have  assessed  and  developed  the  use  of  geostatistical  methods  to  quantify  the  heterogeneity  and 
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structure  seen  in  the  GPR  image;  our  hypothesis  being  that  this  gives  us  a  valid  description  of  the 
heterogeneity  of  hydraulic  properties  of  the  subsurface.  We  analyze  the  GPR  data  by  constructing 
an  experimental  variogram,  where  our  parameter  of  interest  is  the  amplitude  of  the  reflectors  in  the 
GPR  data,  and  then  fit  the  experimental  variograms  with  standard  variogram  models. 

In  order  to  determine  the  link  between  the  spatial  variability  seen  in  the  GPR  image  and  the  tree 
spatial  variability  in  hydraulic  properties  of  the  subsurface,  we  conducted  a  cliff-face  experiment. 
In  this  experiment  we  collected  GPR  data  across  the  top  of  the  cliff  and  could  examine  and 
photograph  the  imaged  sedimentary  sequence,  which  was  a  package  of  alternating  coarse-grained 
sand  and  fine-grained  silt.  We  used  a  geostatistical  analysis  of  the  photograph  of  the  cliff  to 
quantify  the  spatial  variability  or  correlation  structure  in  grain  size  in  the  sediments.  The 
correlation  structure  we  obtained  from  the  GPR  image  was  in  excellent  agreement  with  that  from 
the  photograph  indicating  that  the  GPR  image  had  captured  information  about  the  spatial 
distribution  of  the  coarse-grained/fine-grained  sequence. 

In  an  extension  of  this  study,  experimental  variograms  have  been  obtained  from  GPR  data  from  a 
number  of  different  sedimentary  environments.  The  experimental  variograms  are  all  of  excellent 
quality,  and  have  led  us  to  conclude  that  the  geostatistical  analysis  of  GPR  data  does  provide  a  way 
of  quantifying  spatial  variability  in  a  wide  range  of  geologic  settings.  In  addition  we  find 
indications  that  the  correlation  structure  in  the  GPR  image  is  related  to  the  processes  in  the 
depositional  environment,  leading  us  to  suggest  that  there  can  be  characteristic  GPR  “signatures” 
associated  with  different  environments. 

We  conclude  that  the  geostatistical  analysis  of  GPR  data  from  a  site  can  be  used  for  site 
characterization  where  information  is  needed  about  the  nature  and  extent  of  the  spatial  heterogeneity 
of  the  subsurface.  We  have  developed  an  interactive  computer  program,  included  in  the  insert  at  the 
back  of  this  report,  that  makes  it  possible  to  select  and  analyze  regions  within  a  GPR  image,  and 


model  the  resulting  variogram.  This  could  be  used  with  GPR  data  from  any  site  as  a  means  of 
planning  for  or  integrating  data  obtained  through  drilling  and  direct  sampling. 

The  second  aspect  of  our  research  has  involved  a  study  of  the  direct  relationship  between  the 
measured  values  of  the  dielectric  constant  of  a  volume  of  the  subsurface  and  hydraulic  properties. 
Estimates  of  dielectric  constant  can  be  extracted  from  GPR  data  and  can  also  be  obtained  through 
the  use  of  time  domain  reflectometry.  These  estimates  of  dielectric  constant  can  then  be  used  to 
estimate  hydraulic  properties  through  the  use  of  relationships  between  dielectric  and  hydraulic 
properties.  While  previous  research  has  addressed  these  relationships  for  homogeneous  materials, 
we  have  extended  this  work  to  consider  the  relationship  between  these  properties  in  heterogeneous 
materials.  We  have  considered  a  specific  form  of  heterogeneity,  layering,  which  is  very  likely  to 
be  encountered  in  measurements  of  the  subsurface. 

Theoretical  modeling  of  the  relationship  between  measured  dielectric  constant  and  water  content 
has  demonstrated  the  errors  that  can  result  if  the  sample  size  and  contained  heterogeneity,  such  as 
layering,  are  not  taken  into  account.  In  modeling  the  measured  dielectric  constant,  the  critical 
parameter  is  the  thickness  of  the  layers  relative  to  the  wavelength  of  the  measurement  In  one  limit, 
where  the  wavelength  is  less  than  the  layer  thickness,  ray  theory  is  a  valid  description  of  the 
electromagnetic  response  of  the  material;  in  the  other  limit,  where  the  wavelength  is  greater  that  the 
layer  thickness,  an  effective  medium  approximation  must  be  used.  The  errors  that  result  from 
simply  assuming  a  homogeneous  system  become  pronounced  due  to  the  saturation  heterogeneity 
that  can  develop. 

In  addition  to  theoretical  and  numerical  modeling  of  electromagnetic  wave  propagation  in  layered 
materials,  laboratory  measurements  of  a  sand/silt  system  have  been  conducted.  In  these 
experimental  results  we  can  see  very  clearly  the  predicted  transition  from  ray  theory  to  effective 
medium  theory.  The  results  show  the  complex  nature  of  the  relationship  between  dielectric 


properties  and  hydraulic  properties  in  heterogeneous  materials,  and  highlight  the  need  to  ensure 
that  there  is  a  valid  theoretical  basis  for  the  relationship  that  is  used  to  extract  hydraulic  information 
from  dielectric  measurements. 

In  the  use  of  dielectric  measurements  for  subsurface  characterization,  the  level  of  accuracy  in 
interpretation  of  the  dielectric  data  can  be  significantly  improved  if  additional  parameters  are 
measured.  A  methodology  has  been  developed  that  allows  us  to  extract  electrical  conductivity  from 
GPR  data;  knowledge  of  electrical  conductivity  can  assist  in  constraining  the  interpretation  of  the 
GPR  data  in  terms  of  hydraulic  properties  of  the  subsurface.  Given  that  this  requires  the  collection 
of  no  additional  data,  we  conclude  that  estimates  of  electrical  conductivity  could  be  easily  obtained 
at  a  site  to  improve  the  use  of  GPR  to  estimate  subsurface  hydraulic  properties. 

This  research  has  significantly  advanced  our  understanding  of  how  best  to  use  dielectric 
measurements,  from  GPR  data,  to  both  estimate  the  magnitude  of  hydraulic  properties  and  to 
characterize  spatial  heterogeneity;  both  important  aspects  of  site  characterization.  While  much  of 
this  research  has  involved  the  study  of  fundamental  scientific  issues,  some  of  the  results  of  the 
research  can  be  readily  applied,  now,  to  site  characterization  problems.  In  particular  we  highly 
recommend  that  geostatistical  analysis  of  GPR  images  be  used  as  a  means  of  obtaining  information 
about  the  nature  and  extent  of  spatial  heterogeneity  in  the  subsurface. 

This  report  describes  all  aspects  of  our  three-year  project  Additional  and  more  detailed  information 
can  be  found  in  the  following  publications  and  presentations  that  resulted  from  this  project: 
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1998. 
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INTRODUCTION 

A  critical  step  in  modeling  groundwater  flow  and/or  contaminant  transport  at  a  site  is  to  obtain  an 
accurate  3-D  model  of  the  hydraulic  conductivity  of  the  subsurface.  Past  practice  has  often 
involved  assigning  average  values  of  hydraulic  conductivity  K  to  simple  layered  models  of  the 
subsurface.  However  it  now  widely  recognized  that  using  averaged  values  of  K  is  insufficient; 
the  spatial  variability  in  K,  the  natural  heterogeneity  of  the  subsurface,  has  a  significant  effect  on 
groundwater  flow  and  is  an  important  parameter  in  modeling  contaminant  transport. 

Of  interest  in  this  study  has  been  the  use  of  measurements  of  dielectric  properties  of  the  subsurface 
to  characterize  the  spatial  variability  in  hydraulic  properties  of  the  subsurface.  Measurements  of  in 
situ  dielectric  properties  can  be  made  using  a  ground  penetrating  radar  (GPR)  system.  GPR 
provides  us  with  an  image  of  the  subsurface,  which  is  actually  an  image  of  the  dielectric  properties 
of  the  subsurface.  The  GPR  image  is  obtained  by  sending  a  pulse  of  electromagnetic  energy  into 
the  ground  and  recording  the  amplitude  of  energy  reflected  back  to  the  surface;  energy  is  reflected 
back  from  interfaces  across  which  there  is  a  change  in  dielectric  properties.  Seen  in  the  image  is 
the  location  of  interfaces,  such  that  the  GPR  image  is  actually  an  image  of  the  dielectric  properties 
of  the  subsurface.  One  focus  of  our  research  has  been  to  work  with  GPR  data  and  to  develop  data 
analysis  methods  that  allow  us  to  extract  information  about  the  spatial  variability  in  hydraulic 
properties.  Specifically  we  have  assessed  die  use  of  geostatistical  methods  to  quantify  the 
heterogeneity  and  structure  seen  in  the  GPR  image;  our  hypothesis  being  that  this  gives  us  a  valid 
description  of  the  heterogeneity  of  hydraulic  properties  of  the  subsurface.  This  research  is 
described  in  section  2  of  this  report. 

The  second  aspect  of  our  research,  described  in  section  3,  has  involved  the  direct  relationship 
between  the  measured  values  of  the  dielectric  constant  of  a  volume  of  the  subsurface  and  hydraulic 


properties.  Estimates  of  dielectric  constant  can  be  extracted  from  GPR  data  and  can  also  be 
obtained  through  the  use  of  time  domain  reflectometry.  These  estimates  of  dielectric  constant  can 
then  be  used  to  estimate  hydraulic  properties  through  the  use  of  relationships  between  dielectric  and 
hydraulic  properties.  While  previous  research  has  addressed  these  relationship  between  dielectric 
properties  and  hydraulic  properties  for  homogeneous  materials,  we  have  extended  this  work  to 
consider  the  relationship  between  these  properties  in  heterogeneous  materials.  We  have  considered 
a  specific  form  of  heterogeneity,  layering,  which  is  very  likely  to  be  encountered  in  measurements 
of  the  subsurface. 

In  the  use  of  dielectric  measurements  for  subsurface  characterization,  the  level  of  accuracy  in 
interpretation  of  the  dielectric  data  can  be  significantly  improved  if  additional  parameters  are 
measured.  In  section  4  we  describe  a  way  in  which  electrical  conductivity  can  be  extracted  from 
GPR  data;  knowledge  of  electrical  conductivity  can  assist  in  constraining  the  interpretation  of  the 
GPR  data  in  terms  of  hydraulic  properties  of  the  subsurface. 

As  our  understanding  of  the  dielectric  properties  of  geologic  materials  improves,  so  too  does  our 
ability  to  use  the  measurement  of  these  properties  as  an  effective  means  of  subsurface 
characterization.  We  have  found  that  measurement  of  dielectric  properties  can  be  used  to  obtain 
information  about  both  the  magnitude  and  spatial  heterogeneity  of  hydraulic  properties.  The 
measurement  of  dielectric  properties,  by  the  collection  of  GPR  data  at  a  site,  can  thus  provide 
valuable  information  allowing  us  to  more  accurately  model  the  complex  scale-dependent  hydraulic 
properties  of  the  subsurface. 
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2 

Geostatistical  Analysis  of  Ground  Penetrating  Radar  Data: 

A  Means  of  Characterizing  the  Heterogeneity  of  the  Subsurface 

Introduction 

A  critical  step  in  modeling  contaminant  transport  or  remediation  at  a  site  is  obtaining  a  realistic  3-D 
model  of  a  region  of  the  subsurface.  Given  the  heterogeneity  of  the  geological  environment,  it  will 
never  be  possible  to  collect  sufficient  data  through  drilling  and  direct  sampling  to  adequately 
characterize  the  spatial  variability  that  exists  in  the  subsurface.  The  question  becomes  how  best  to 
"fill  in"  the  region  between  data  points,  or  describe  the  region  at  a  scale  less  than  the  scale  of  our 
sampling.  A  common  approach  is  to  use  geostatistical  methods. 

Geostatistics  provides  a  mathematical  framework  that  can  be  used  to  describe  how  some  property 
or  parameter  varies  spatially;  these  spatial  statistics  can  then  be  used  in  the  development  of  the 
model  of  the  subsurface.  Some  information  about  spatial  statistics  can  be  obtained  from  borehole 
data  or  cores,  but  due  to  the  relatively  large  distances  between  boreholes,  this  usually  only 
provides  a  good  measure  of  variation  in  the  vertical  direction.  To  obtain  a  description  of  lateral 
variation  data  can  be  collected  from  outcrops,  selected  as  representative  of  the  subsurface  geology, 
but  this  is  extremely  time  consuming,  and  usually  limited  to  the  two  dimensions  of  the  exposed 
section. 

We  suggest  that  geostatistical  analysis  of  GPR  data  can  be  used  to  describe  the  spatial  variation  of 
the  subsurface.  Our  working  hypothesis  is  that  the  links  between  sedimentary  and  hydrogeologic 
properties,  and  between  sedimentary  and  dielectric  properties,  can  be  exploited  to  infer  an 


underlying  relationship  between  the  spatial  distribution  of  the  radar  reflectors  seen  in  a  GPR 
section  and  the  spatial  distribution  of  hydrogeologic  properties. 


In  our  work  we  use  the  semivariogram,  a  plot  which  illustrates  the  way  in  which  the  difference 
between  data  values  is  related  to  the  distance  between  the  data  values.  The  experimental 
semivariogram  is  described  by  the  following  equation  (Deutsch  and  Joumel,  1992): 


Y(h)  = 


1 

2N(h) 


N(h)  2 

X[z(xi)~z(xi  +h)] 


i*=l 


(1) 


where  h  is  the  lag,  or  separation  vector,  between  two  data  points,  z(x)  and  z(x+h)  and  N  is  the 
number  of  data  pairs  used  in  each  summation.  In  general,  the  semivariogram  function  y  increases 

with  increasing  lag,  leveling  off  at  the  "sill"  when  the  distance  between  the  data  points  becomes  so 
large  that  the  data  values  are  no  longer  correlated.  The  value  of  the  lag  at  this  point  is  referred  to  as 
the  range  or  the  correlation  length  of  the  data  set.  In  most  cases  a  semivariogram  is  directional,  that 
is,  the  vector  h  lies  in  a  specified  direction.  This  direction  can  be  set  to  any  direction  of  interest, 
but  is  generally  that  of  maximum  correlation.  As  our  objective  is  to  determine  the  correlation 
length  of  the  radar  reflectors,  we  use  the  amplitude  values  recorded  in  the  radar  traces  as  our 
parameter  z.  We  obtain  the  experimental  variogram  for  the  GPR  image  and  suggest  that  this  is 
representative  of  the  correlation  structure  of  the  subsurface;  i.e.  the  GPR  image  provides  a  way  of 
quantifying  the  nature  of  the  spatial  heterogeneity  in  the  subsurface  properties. 

Classification  of  Radar  Units 

An  underlying  assumption  in  the  calculation  of  semivariograms  is  that  the  data  conform  to  a 
requirement  called  stationarity.  Stationarity  means  that  any  subset  of  the  data  will  have  the  same 


10 


statistical  description  as  any  other  subset  This  is  of  concern  in  most  geostatistical  problems  as  the 
earth  is  not  inherently  stationary,  with  different  geologic  units  displaying  very  different  spatial 
variability.  For  the  purposes  of  our  study,  this  implies  that  we  must  separate,  and  analyze 
individually,  data  obtained  from  different  sedimentary  units. 

Working  with  a  large  data  set  of  GPR  images,  collected  in  the  area  of  the  Brookswood  aquifer, 
southwestern  British  Columbia,  we  developed  a  methodology  for  defining  “radar  elements”;  this 
is  described  in  detail  in  Rea’s  PhX).  thesis  (appendix  A)  and  Rea  and  Knight  (1995).  The  method 
involves  first  developing  a  classification  scheme  that  makes  it  possible  to  relate  different  regions  in 
the  GPR  images  to  the  different  sedimentary  or  hydrogeologic  units.  This  was  done  primarily 
using  the  geometry  of  packages  of  reflectors  (e.g.  “channel-shaped)  and  the  length  and  dip  of  the 
internal  reflectors  (e.g.  “long,  subhorizontal  reflectors”).  Well  data  allowed  us  to  relate  these  GPR 
characteristics  to  the  different  hydrogeologic  units. 

Method  of  Geostatistical  Analysis 

In  our  geostatistical  analysis  of  the  GPR  data,  our  objective  is  to  determine  the  correlation  length  of 
the  radar  reflectors.  The  data  values  we  use  to  construct  the  experimental  semivariogram  are  the 
amplitude  values  recorded  in  the  radar  traces.  We  first  correct  for  a  violation  of  the  stationarity 
assumption  introduced  by  the  GPR  method  due  to  the  decay  of  amplitude  down  a  radar  trace. 
Uncorrected,  the  data  at  early  times,  having  much  higher  amplitudes,  would  dominate,  the 
semivariogram  calculation.  We  first  correct  the  data  for  any  time  shifts  due  to  elevation  changes  or 
problems  with  the  electronics,  then  use  an  automatic  gain  control  with  a  short  window  length  to 
make  the  amplitudes  of  all  the  reflectors  approximately  equal.  The  high  amplitude  air  and  ground 
arrivals  at  the  start  of  each  trace  are  removed  by  simply  cropping  the  data.  A  further  processing 
step,  migration  of  the  data,  might  be  needed  to  recover  the  true  correlation  length  in  areas  with 
steeply  dipping  reflectors;  migration  was  found  to  have  little  effect  on  the  analysis  of  the  data  from 


our  study  areas.  As  a  final  step  before  calculating  a  semivariogram,  the  radar  Hat  a  must  be 
converted  from  a  time  section  to  a  depth  section  so  that  the  lag  is  a  physical  length  scale.  This  is 
done  using  the  velocity  calculated  from  a  common  midpoint  radar  survey  conducted  at  the  site 
under  investigation. 

In  our  study  of  the  geostatistical  analysis  of  GPR  data,  we  began  by  using  the  geostatistical 
software  package  GSLIB  (Deutsch  and  Joumel,  1992).  This  was  however  quite  a  tedious 
process.  Several  different  FORTRAN  routines  from  GSLIB,  in  a  non-graphical  environment,  had 
to  be  used  and  it  was  only  possible  to  compute  the  semivariogram  in  one  direction.  As  a  result  we 
developed  a  program  here  referred  to  as  GSWIN,  a  Windows  95  based  program  for  performing 
geostatistical  analysis  on  ground  penetrating  radar  data;  a  copy  is  included  on  the  disc  submitted 

with  this  report.  GSWIN  is  capable  of  displaying  the  variation  in  the  semivariogram  function  y  in 
all  directions  and  does  so  employing  a  graphical  user  interface. 


GSWIN  can  import  three  types  of  data  file  formats: 

1.  GPR  data  collected  and  processed  using  Sensors  and  Software’s  PULSE  EKKO  GPR  package. 
The  resulting  text  file  can  then  be  loaded  into  GSWIN  (Figure  2.1).  A  wave  velocity  must  be 
specified  by  the  user  in  order  to  determine  the  sampling  rate  in  the  y-direction. 

2.  A  SCION  Image  (SCION  Corporation,  1997)  picture  which  has  been  saved  in  the  text  format 
The  sampling  rates  in  the  x  and  y  directions  must  be  entered  by  the  user  in  units  of  metres/pixel. 

3.  An  SGSIMM  simulation.  SGSIMM  takes  an  input  semivariogram  and  simulates  a  data  set. 
This  was  primarily  used  for  testing  whether  the  program  was  working  properly  but  remains  part  of 
the  program  as  an  option.  A  matrix  size  must  be  entered  by  the  user  to  specify  the  size  of  the  data 
set  because  SGSIMM  outputs  numbers  in  one  column  only. 

Once  the  data  are  in  the  main  window,  a  variogram  surface  can  be  computed.  This  gives  the  spatial 
variation  of  the  data  in  all  directions;  i.e.  the  lag  vector  h  can  have  any  direction.  This  computation 


Figure  2.1.  Data  imported  from  PULSE  EKKO  format. 


is  done  by  comparing  all  points  in  the  grid  to  each  other.  The  user  can  choose  to  select  a 
rectangular  or  irregular  region  on  which  to  perform  the  analysis  if  the  entire  data  set  is  not  desired. 
In  this  way  it  is  possible  to  separate  for  analysis  regions  in  the  GPR  image  that  represent  different 
geologic  units.  These  options  appear  in  the  toolbar  as  well  the  pull-down  menus.  Once  the  mode 
is  set,  the  user  draws  the  region  by  holding  the  left  mouse  button  and  moving  the  mouse.  When  a 

region  has  been  selected,  clicking  the  semivariogram  toolbar  button  (y)  gives  the  user  a  list  of 

options  (Figure  2.2)  for  computing  the  variogram  surface.  The  user  can  either  change  the  number 
of  lags,  which  determines  the  number  of  data  points  in  the  plot,  the  lag  spacing,  which  indicates 
how  many  metres  each  pixel  in  the  semivariogram  represents,  or  leave  them  on  the  defaults.  These 
defaults  are  the  x  and  y  sampling  rates  for  the  data  set,  and  30  pixels  in  each  direction.  Once  these 
are  specified  the  program  computes  a  two-dimensional  semivariogram  on  the  data  (Figure  2.3). 
This  window  is  split  in  two,  one  containing  the  actual  semivariogram,  the  other  showing  the 
number  of  pixels  and  lag  spacing.  In  order  to  calculate  the  one  dimensional  semivariogram  in  a 
particular  direction,  which  is  the  ultimate  goal,  the  user  clicks  on  a  portion  of  the  variogram  surface 
with  the  left  mouse  button  and  a  line  is  drawn  from  the  center  out  to  that  point  and  the  angle  is 
calculated.  This  line  can  then  be  moved  by  moving  the  mouse  and  holding  down  the  left  button. 
When  the  angle  is  to  the  user’s  satisfaction,  the  1-D  semivariogram  is  calculated  (Figure  2.4)  when 
the  user  releases  the  left  mouse  button. 

An  experimental  semivariogram  results  in  values  for  y  only  at  certain  lags,  which  are  determined  by 
the  sampling  interval  of  the  data.  Semivariogram  models  are  used  to  provide  an  analytic 
description  of  the  experimental  semivariogram  so  that  y  can  be  calculated  for  any  lag  h.  This  is 

needed  in  order  to  use  the  semivariogram  for  simulation  algorithms  and  to  provide  a  better  ftsthmtp. 
of  the  true  correlation  length  of  the  data  set.  In  Figure  2.4  the  1-D  semivariogram  window  is  a 
split  window  showing  the  semivariogram  in  the  left  half  and  the  fitting  parameters  in  the  right. 


Figure  2.2.  The  options  for  the  variogram  computation. 


Figure  2.3.  Two  dimensional  variogram.  The  angle  of  the  dark  ellipse  indicates  the  maximum 
correlation  direction. 


Figure  2.4.  A  one  dimensional  variogram.  The  exponential  model  was  the  best  fit  in  this  case. 


Here  the  user  can  determine  which  model  best  fits  the  data.  The  models  which  the  program 
supports  are  as  follows  (Jian  et  al,  1996): 


1. Spherical 


2.Exponential 


3.Power 
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C0  =  nugget 
a  =  range 
C  =  sill 


In  order  to  determine  the  best-fit  model,  the  user  can  press  the  ‘best  fit’  button  shown  in  Figure  2.4 
and  the  program  will  perform  a  nonlinear  least-squares  fit  using  each  model.  The  one  with  the 

lowest  resulting  %2  is  chosen  as  the  best-fit.  The  fitted  semivariogram  is  drawn  in  the  left  window 

using  a  green  line.  A  nested  model  (linear  combination)  of  up  to  three  different  models  may  also 
be  used  by  specifying  the  models  in  the  boxes  and  choosing  ‘fit’.  Hie  initial  guesses,  which 
appear  in  the  top  half  of  the  window,  can  be  changed  by  the  user,  though  the  defaults  are  often 


19 


close  to  the  final  value  anyway.  The  defaults  are  chosen  by  an  algorithm  which  guesses  at  the 
range  by  estimating  where  the  graph  begins  to  level  out.  The  sill  is  chosen  by  averaging  all  points 
after  this  estimated  range.  The  user  can  also  select  an  individual  model  and  press  the  ‘fit’  button  to 
determine  the  parameters  for  a  particular  model.  Once  the  analysis  is  complete,  the  user  can  choose 
to  save  the  semivariogram  along  with  all  the  fitted  parameters  so  it  may  be  loaded  into  a  graphing 
package,  such  as  Excel.  The  user  can  also  reduce  the  number  of  points  to  be  fit  by  clicking  the  left 
mouse  button  on  the  variogram  window  and  moving  the  cutoff  line  to  the  desired  position. 

This  program  was  written  using  Microsoft  Visual  C++.  The  nonlinear  least  squares  fitting  routine 
uses  the  Levenberg-Marquardt  technique  as  implemented  in  Numerical  Recipes  (Press  et  al,  1992). 
The  calculation  of  the  2-D  semivariogram  was  based  on  a  routine  from  GSLIB  which  was 
converted  to  C++  from  its  original  FORTRAN. 

Calibration  Field  Experiment 

Our  working  hypothesis  is  that  the  correlation  structure  determined  for  the  GPR  image  is  a  valid 
description  of  the  correlation  structure  of  the  hydraulic  properties  of  the  subsurface.  This  assumes 
that  the  spatial  variation  in  dielectric  properties,  seen  in  the  GPR  data,  is  an  accurate  image  of  the 
spatial  variation  in  hydraulic  properties  of  the  subsurface.  The  critical  issue  which  must  be 
addressed  then  is  the  relationship  between  the  GPR  image  and  the  properties  of  the  subsurface. 
Determining  this  link,  between  the  GPR  image  and  the  "geological  reality"  is  a  key  part  of  our 
research.  The  approach  that  we  are  taking  is  to  conduct  cliff  face  studies:  GPR  data  are  collected 
along  the  top  of  the  cliff  and  the  imaged  section  can  be  seen  and  sampled.  This  offers  an 
outstanding  opportunity  to  determine  what  is  really  imaged  in  the  GPR  section;  i.e.  to  "ground 
truth"  the  GPR  data. 


In  one  cliff  face  experiment  (described  in  Rea  and  Knight,  1998)  we  compared  the  geostatistical 
analysis  of  a  photographic  image  and  a  GPR  image  of  a  sequence  of  alternating  coarse  sand,  and 
fine  sand  and  silt.  The  digital  photograph  of  the  face  captured  information  about  the  spatial 
distribution  of  coarse  grained  and  fine  grained  beds  on  the  basis  of  gray  scale.  There  was  excellent 
agreement  between  the  variograms  from  the  photograph  and  the  radar  data,  in  determining  both  the 
maximum  correlation  direction  and  the  correlation  length.  These  results  led  us  to  conclude  that  die 
GPR  data  did  in  fact  image  the  spatial  distribution  of  these  two  lithologies  and  could  be  used  to 
quantify  both  the  correlation  direction  and  length  in  this  sedimentary  unit  In  this  example  die 
spatial  variation  in  dielectric  properties  in  the  subsurface  was  closely  related  to  the  spatial  variation 
in  grain  size;  this  is  very  reasonable  given  what  is  known  about  the  relationship  between  dielectric 
properties  and  sedimentary  properties. 


Analysis  of  GPR  Data  from  Selected  Depositional  Environments 

In  collaboration  with  Harry  Jol,  University  of  Wisconsin  Eau-Claire,  we  conducted  geostatistical 
analyses  of  GPR  images  from  a  variety  of  selected  depositional  environments.  We  had  two 
principal  objectives  in  this  work.  The  first  was  to  determine  whether  the  geostatistical  analysis  of 
GPR  data,  that  we  had  shown  to  be  successful  in  the  cliff-face  example  described  in  Rea  and 
Knight  (1998),  would  be  successful  in  a  wide  range  of  sedimentary  environments.  Our  second 
objective  was  to  investigate  a  more  fundamental  issue:  We  suggest  that  the  physical  and  chemical 
processes  acting  in  any  given  depositional  environment  should  be  responsible  for  determining  the 
spatial  heterogeneity  in  material  (e.g.  hydraulic)  properties  in  that  environment  We  want  to  test 
the  idea  that  this  will  lead  to  different  depositional  environments  having  different  but 
“characteristic”  geostatistical  signatures.  That  is,  the  semivariograms  that  represent  sedimentary 
packages  from  braided  streams  will  all  have  certain  features  in  common,  and  will  differ  in 
predictable  ways  from  semivariograms  for  sedimentary  packages  from  coastal  spits.  The 
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characteristic  features  might  be  the  correlation  length  or  the  functional  form  of  the  semivariogram 
(exponential,  spherical,  nested  etc.). 

The  depositional  environments  that  we  have  examined  are  listed  in  Table  2.1.  In  our  work  we 
have  achieved  our  first  objective  -  the  geostatistical  analysis  of  GPR  data  can  be  conducted  in  a 
wide  range  of  depositional  environments.  In  all  cases  we  have  obtained  outstanding  experimental 
semivariograms  that  are  well-modeled  by  the  standard  semivariograms  models.  With  respect  to 
our  second  objective,  we  do  observe  differences  in  the  geostatistical  “signature”  of  different 
environments,  and  can  make  some  general  conclusions  about  the  relationship  between  depositional 
processes  and  the  resulting  spatial  statistics.  However,  further  work  is  required  to  develop  a  more 
fundamental  understanding  of  the  link  between  depositional  environments  and  their  geostatistical 
character.  The  preliminary  results  from  this  work  are  described  in  detail  in  Knight  et  al.  (1997); 
the  final  results  were  presented  recendy  at  the  Hazardous  Substances  Research  Conference  (Knight 
et  al.,  1998)  and  are  in  preparation  as  a  paper  to  be  submitted  to  the  journal  Geophysics.  Our 
principal  findings  are  briefly  described  below. 

Table  2.2  lists  the  results  of  the  geostatistical  analysis  of  the  GPR  images.  We  find  correlation 
lengths  ranging  from  3m  to  53  m  and  find  that  an  exponential  model  best  describes  the  majority  of 
the  data.  As  one  example  of  the  differences  we  see  in  the  geostatistical  signature  of  different 
environments,  we  can  compare  the  GPR  image  of  a  fan-foreset  delta  and  the  GPR  image  of  a  braid 
delta.  As  can  be  seen  in  the  GPR  image  of  the  fan-foreset  delta  (Figure  2.5)  this  environment 
contains  long  continuous  dipping  features  that  are  found  to  have  a  maximum  correlation  direction 
of  20°  to  the  northwest..  These  features  in  the  GPR  data  are  interpreted  to  represent  the  foreset 
beds  within  the  delta.  In  Figure  2.6  is  shown  the  experimental  semivariogram  in  that  direction, 
with  a  correlation  length  of  53  m.  In  contrast  to  the  depositional  environment  for  a  fan-foreset 
delta,  which  tends  to  produce  and  preserve  long,  continuous  bedding  planes,  a  braid  delta  is 
formed  in  an  environment  dominated  by  high  fluctuations  in  energy  level.  As  can  be  seen  in  the 


Table  2.1.  Depositional  Environment 
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Table  2.2.  Semivariogram  analysis:  parameters  and  results 


Environment  & 
Location 

Facies 

Angle 

_0 

model 

range 

(m) 

DELTAS 

b 

i  U 

is 

c 

Box  Elder 
Creek 

Foreset 

20  NW 

exp. 

gauss. 

53 

3 

BRAID 

American 
Forks  R. 

Cut  and  Fill 

0.5 

exp. 

3 

Athabasca 

River 

A  (cut  and  fill) 

1  W-NW 

exp. 

6 

B  (cut  and  fill) 

1  E-SE 

exp. 

4 

C  (cut  and  fill) 

1  E-SE 

exp. 

5 

D  (cut  and  fill) 

0 

exp. 

4 

WAVE 

William 

River 

Uppershore/ 

Beachface 

IE 

exp. 

7 

Lower-mid  shore 

2E 

pent. 

14 

SPITS 

MARINE 

Willapa 

Bay 

Dip 

Beachface 

1.5  W 

exp. 

24 

Strike 

Beachface 

0 

sph. 

gauss. 

6 

46 

LACUSTRINE 

Beachface  || 

1  SE 

pent. 

18 

Beachface  ± 

0 

sph. 

7 

1 

Stormzone  || 

2  SE 

exp. 

7 

| 

Stormzone  ± 

3  SW 

exp. 

10 

Lakebed  || 

0 

exp. 

43 

Lakebed  J_ 

2  SW 

exp. 

9 

exp.  =  exponential,  gauss.  =  gaussian,  pent.  =  pentaspherical,  sph.=  spherical 
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Box  Elder  Creek  fan-foreset  delta 
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Figure  2.5.  GPR  profile  collected  over  the  fan-foreset  delta  (modified  from  Smith  and  Jol, 
1992a). 
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Figure  2.6  Semivariogram  of  the  GPR  data  collected  over  the  Box  Elder  Creek  delta  foreset 
beds  shown  in  Figure  2.5.  The  data  points  are  the  circles;  the  model  of  the  semivariogram  is  the 
solid  line.  These  data  are  modeled  using  a  nested  exponential  and  gaussian  model  with  ranges  of 
53  m  and  3  m.  This  very  long  range,  53  m,  is  a  measure  of  the  distance  over  which  the  inclined 
strata  are  continuous. 


GPR  image  from  a  braid  delta  (Figure  2.7),  the  section  is  made  up  of  numerous  discontinuous 
features  that  are  interpreted  as  “cut  and  fill”  structures  which  reflect  the  continuous  shifting  of  the 
channel  position.  The  semivariogram  from  this  environment,  shown  in  Figure  2.8,  has  a 
correlation  length  of  3m.  Clearly  we  see  in  the  geostatistical  “signature”  of  these  two  deltas  what 
we  would  expect  given  the  nature  of  the  depositional  environments  -  long  correlation  lengths  in  the 
fan-foreset  delta,  and  much  shorter  correlation  lengths  in  the  braid  delta. 


As  a  second  example  of  differences  in  correlation  structure  due  to  differences  in  depositional 
processes  we  show  in  Figure  2.9  the  GPR  data  from  a  lacustrine  spit.  In  this  case  there  is  a 
change  in  depositional  processes  with  time  -  which  corresponds  to  progradation  of  the  spit  into  the 
lake.  In  analyzing  this  data  set  GSWIN  proved  to  be  extremely  useful,  as  we  could  select  and 
analyze  separately  three  regions  in  the  data,  as  shown  on  the  GPR  section.  Again,  in  this  example 
we  see  differences  in  the  semivariograms  (Figures  2.10  a,b,c)  that  can  be  qualitatively  related  to 
die  depositional  processes.  The  lake  bed  represented  by  long,  horizontal  reflections  has  a 
maximum  correlation  length  of  43m.  The  storm  generated  deposits  reflect  a  high  energy 
environment  and  are  characterized  by  shorter,  dipping  reflectors  with  a  correlation  of  7m.  The 
beach  face  which  corresponds  to  the  zone  of  wave  swash  is  represented  by  long  horizontal 
reflections  with  a  correlation  length  of  18m. 

As  can  be  seen  in  the  above  examples,  the  quality  of  the  all  the  variograms  are  outstanding.  We 
conclude  that  geostatistical  analysis  of  GPR  data  can  be  applied  in  a  variety  of  geological 
environments.  The  examples  also  illustrate  the  concept  which  we  have  begun  to  investigate  -  can 
the  variograms  from  different  environment  be  explained  in  terms  of  the  processes  occurring  in  the 
environments?;  and  is  there  a  characteristic  geostatistical  "signature"  associated  with  specific 
environments?  Such  signatures  would  be  extremely  useful  in  attempting  to  interpret  and  integrate 
data  at  a  site  if  there  was  some  prior  knowledge  of  the  depositional  environment 


American  Fork  River  braid  delta 


Figure  2.7.  GPR  profile  collected  over  the  American  Fork  River  braid  delta  (modified  from  Jol 
and  Smith,  1992). 
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Figure  2.8.  Semivariogram  analysis  of  the  GPR  image  shown  in  Figure  2.7.  The  data  points  are 
the  circles;  the  model  of  the  semivariogram  is  the  solid  line.  These  data  are  modeled  using  an 
exponential  model  with  a  range  of  3  m.  These  short  discontinuous  reflectors  have  been 
interpreted  as  cut  and  fill  structures. 


Sandy  Point  spit:  parallel  to  spit  axis 


Figure  2.9.  GPR  profile  collected  parallel  to  the  spit  axis  over  the  active  Sandy  Point  spit 
(modified  from  Smith  and  Jol,  1992b).  Areas  outlined  indicate  data  used  in  geostatistical 
analysis.  Area  “A”  represents  beachface/upper  shoreface,  area  “B”  represents  storm  influenced 
zone,  and  area  “C”  represents  the  lake  bed. 


Depth  (m) 


30 


Sandy  Point  spit:  parallel  to  spit  axis 
Area  A:  beachface/upper  shoreface 
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Figure  2.10a.  Semivariogram  analysis  of  area  “A”  in  Figure  2.9.  These  data  are  modeled  using 
a  pentaspherical  model  with  a  range  of  1 8  m. 
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Sandy  Point  spit:  parallel  to  spit  axis 
Area  B:  storm  zone 
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Figure  2.10b.  Semivariogram  analysis  of  area  “B”  in  Figure  2.9.  These  data  are  modeled  using 
an  exponential  model  with  a  range  of  7  m. 
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Figure  2.10c.  Semivariogram  analysis  of  area  “C”  in  Figure  2.9.  These  data  are  modeled  using 
an  exponential  model  with  a  range  of  43  m. 
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Conclusions 

The  results  presented  in  this  section  illustrate  the  potential  value  of  using  GPR  images  of  the 
subsurface  to  describe  and  quantify  the  nature  of  the  spatial  heterogeneity  in  the  subsurface.  Such 
information  could  be  used  at  the  start  of  a  site  characterization  exercise  in  order  to  plan  the  spacing 
and  density  required  for  drilling  and  direct  sampling.  In  addition,  the  correlation  structure  imaged 
with  GPR  could  be  used  to  assist  in  the  integration  and  interpretation  of  data  from  a  site. 
Inevitably,  in  the  assessment  of  any  site  there  will  be  a  shortage  of  data.  Information  about  spatial 
variability,  as  can  be  provided  by  GPR  images,  can  assist  in  developing  models  of  the  subsurface 
that  adequately  account  for  the  heterogeneity  in  hydraulic  properties. 
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3 

The  Relationship  Between  Dielectric  Constant 
and  Hydraulic  Properties 

Introduction 

Of  interest  in  our  research  is  an  assessment  of  the  ways  in  which  measurements  of  dielectric 
properties  of  the  subsurface  can  be  used  to  obtain  information  about  hydrogeologic  properties.  In 
the  previous  section  we  described  a  methodology  for  using  GPR  images  (i.e.  dielectric  images)  to 
characterize  the  spatial  heterogeneity  in  subsurface  materials.  In  this  chapter  we  describe  research 
focused  on  using  measurements  of  the  dielectric  constant  of  a  volume  of  the  subsurface  to  obtain 
estimates  of  the  magnitude  of  hydrogeologic  properties.  That  is,  can  we  measure  the  dielectric 
constant  of  an  area  and  from  that  obtain  accurate  estimates  of  water  content,  porosity,  permeability? 

There  are  two  methods  commonly  used  to  determine  the  dielectric  constant  of  a  region  in  tire 
subsurface.  In  time  domain  reflectometry  (TDR)  the  dielectric  constant  is  calculated  from  the  time 
it  takes  an  electromagnetic  (EM)  wave  to  travel  along  the  TDR  probe,  a  wave  guide,  inserted  into 
the  ground.  The  sampled  material  over  which  the  dielectric  constant  is  determined  is  the  volume  of 
material  immediately  adjacent  to  and  along  the  full  length  of  the  TDR  probe.  A  TDR  study  is  most 
commonly  conducted  to  assess  the  variation  in  water  content  in  the  top  few  meters  of  the  earth  and 
typically  covers  an  area  of  tens  of  square  meters  with  sample  spacing  on  the  order  of  a  meter. 

Dielectric  information  can  also  be  obtained  using  a  ground  penetrating  radar  (GPR)  system.  In 
using  GPR,  two  antennae  are  moved  across  the  earth’s  surface:  one  to  transmit  electromagnetic 
energy  and  the  other  to  receive  energy  that  has  been  reflected  back  to  the  surface  from  interfaces 
across  which  there  have  been  changes  in  dielectric  properties.  The  volume  of  the  subsurface 
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sampled  in  GPR  studies  usually  ranges  from  hundreds  to  thousands  of  cubic  meters.  The 
dielectric  constant  of  a  region  determines  the  velocity  of  the  electromagnetic  wave  such  that  GPR 
velocity  data  can  be  used  to  obtain  a  model  of  the  dielectric  constant  of  the  subsurface. 

The  relationship  between  the  measured  dielectric  constant  and  die  hydrogeologic  properties  (water 
saturation,  porosity,  permeability  or  hydraulic  conductivity)  is  die  key  to  our  ability  to  transform 
measured  values  of  dielectric  constant  in  the  subsurface  to  the  hydrogeological  properties  of 
interest  Our  concern  is  that  all  laboratory-derived  relationships  to  date  are  for  homogeneous 
samples.  Similarly,  all  the  theoretical  and  empirical  relationships  commonly  used  to  relate 
dielectric  measurements  to  other  properties  assume  a  homogeneous  material.  The  unique  aspect  of 
our  study  is  an  assessment  of  the  way  in  which  the  presence  of  heterogeneity  changes  the 
relationship  between  dielectric  and  hydrogeologic  properties. 

While  we  plan  eventually  to  assess  the  relationships  between  dielectric  constant  and  various 
hydrogeologic  properties,  we  selected  for  this  study  the  dielectric  constant  -  water  content 
relationship,  as  measurements  of  dielectric  constant  are  most  commonly  conducted  to  monitor  or 
estimate  water  content.  In  addition,  we  elected  to  start  with  a  specific  form  of  heterogeneity  - 
layering.  This  is  a  form  of  heterogeneity  commonly  encountered  as  most  TDR  measurements  and 
many  GPR  measurements  are  made  in  sedimentary  environments  where  the  soils,  sediments,  and 
rocks  arc  complex  interlayering  of  different  lithologies.  These  layers  commonly  parallel  the 
surface  of  the  earth  and  can  range  in  thickness  from  millimeters  to  many  meters. 

We  divided  this  study  into  three  parts  -  described  in  the  three  sections  below.  First  the  theoretical 
differences  between  the  dielectric  constant-water  content  and  dielectric  constant-water  saturation 
relationships  for  homogeneous  systems  and  layered  systems  were  explored.  We  then  numerically 
modeled  the  propagation  of  electromagnetic  waves  through  layers  of  contrasting  dielectric 
constants.  Finally,  laboratory  measurements  were  made  of  electromagnetic  velocity  through 


36 


layers.  From  these  three  aspects  of  our  study  we  have  developed  an  improved  understanding  of 
TDR  and  GPR  measurements  in  layered  materials  and  suggest  how  this  can  lead  to  improved 
accuracy  in  the  use  of  dielectric  data  to  estimate  hydrogeologic  parameters. 

Theoretical  Study  of  Properties  of  Layered  Materials 

Introduction 

This  theoretical  study  examined  the  relationships  between  dielectric  constant  and  water  content  and 
water  saturation  in  layered  materials.  Our  main  objective  in  this  study  was  to  determine  the  level  of 
inaccuracy  in  estimates  of  water  content  or  saturation  when  the  sampled  system  is  incorrecdy 
assumed  to  be  homogeneous.  The  assumption  of  homogeneity,  when  layering  is  present,  can 
introduce  significant  error  into  the  determined  water  content  or  saturation  levels.  A  complete 
description  of  this  study  is  contained  in  Chan  and  Knight  (1998)  recently  submitted  to  Water 
Resources  Research;  the  main  results  are  summarized  below. 

Description  of  Models  and  Theories 

Three  model  systems  composed  of  sand  and  clay,  sand  and  silty  clay,  and  sand  and  silt  loam  were 
considered.  In  each  of  these  systems  the  composition  of  the  solid  phase  and  the  water/air  content 
were  varied.  We  first  describe  each  of  these  systems  as  a  homogeneous  mixture  of  components 
and  then  as  a  set  of  layers.  For  each  case,  for  each  system,  we  theoretically  predict  the  relationship 
between  dielectric  constant  and  water  content  and  dielectric  constant  and  water  saturation. 

In  theoretically  modeling  the  homogeneous  system,  a  specific  form  of  binary  mixtures,  previously 
studied  by  Nur  et  al.  (1991)  and  Knoll  and  Knight  (1996)  was  modeled.  Water  and  air  were 
assumed  to  be  distributed  evenly  such  that  the  average  or  global  water  content  or  saturation  of  the 
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sand-fine  mixture  can  be  considered  to  be  uniform  throughout  the  system.  The  relationships  used 
to  relate  dielectric  constant  to  water  content  and  water  saturation  were  the  TP  (time  propagation) 
model  and  the  Topp  equation.  These  two  models,  which  are  commonly  used  to  interpret  measured 
dielectric  constants,  assume  a  homogeneous  mixture  of  air,  water,  and  solid. 

The  layered  systems  which  were  theoretically  modeled  were  composed  of  the  sand  and  fine 
fraction  used  in  the  homogeneous  mixtures  but  arranged  in  distinct  layers.  Each  layer  was 
composed  of  a  single  sediment  type.  The  volume  fraction  of  fines  was  varied  from  0  to  1  and  the 
global  water  content  from  0  to  the  value  obtained  when  the  pore  space  of  the  entire  system  is  filled 
with  water.  A  critical  issue  in  the  layered  model  was  the  heterogeneity  that  will  exist  in  the 
distribution  of  water  in  the  system.  Under  natural  conditions  differences  in  water  content  exist 
between  sand  and  fine  layers  which  in  turn  translate  into  differences  in  water  saturation  for  the 
layers.  In  order  to  look  at  the  dependence  of  dielectric  constant  on  global  water  content  or 
saturation  in  a  layered  sand-fine  system,  the  water  content  and  saturation  of  the  individual  layers 
were  determined  using  Clapp  and  Homberger’s  (1978)  capillary  pressure-saturation  relationship. 
The  dielectric  constant  of  the  individual  layers  were  predicted  using  either  the  TP  model  or  the 
Topp  equation.  These  values  were  then  used  in  calculating  the  average  dielectric  constant  and  the 
global  water  content  and  water  saturation  for  the  entire  system. 

The  orientation  of  the  layering  and  the  thickness  of  the  layers  were  additional  variables  in  the 
modeling.  We  consider  an  electromagnetic  (EM)  wave  traveling  through  the  layered  material,  both 
parallel  to  and  perpendicular  to  the  layering.  We  consider  two  end-member  cases  where  the 
wavelength  of  the  EM  wave  is  much  greater  than,  and  much  less  than,  the  thickness  of  the  layers. 
In  the  first  case,  the  average  dielectric  properties  are  modeled  using  an  effective  medium  theory;  in 
the  second  case,  ray  theory  is  used. 


Data/Results  and  Discussion 


Plots  of  all  the  theoretical  modeling  results  are  found  in  Chan  and  Knight  (1998).  Figure  3.1 
shows  an  example  of  the  comparison  of  results  from  the  theoretical  calculations:  the  sand-clay 
system  is  shown  at  a  fine  fraction  of  0.5.  Both  the  dielectric  constant-water  content  relationship 
(upper  plot)  and  the  dielectric  constant-water  saturation  relationship  (lower  plot)  are  considered. 
The  relationships  for  homogeneous  systems  are  theoretically  predicted  using  both  TP  and  Topp. 
For  the  layered  systems,  the  dielectric  constants  of  individual  layers  were  calculated  using  the  TP 
model.  We  then  model  three  cases  of  layered  systems: 

1)  wave  propagation  direction  is  perpendicular  to  the  layering,  and  the  layer  thickness  is  much  less 
than  the  wavelength;  this  is  referred  to  as  “perpendicular  EMT”. 

2)  wave  propagation  direction  is  perpendicular  to  the  layering,  and  the  layer  thickness  is  much 
greater  than  the  wavelength;  this  is  referred  to  as  “perpendicular  ray”. 

3)  wave  propagation  direction  is  parallel  to  the  layering,  and  the  layer  thickness  is  much  less  than 
the  wavelength;  this  is  referred  to  as  “parallel  EMT’. 

The  main  objective  in  this  study  was  to  gain  some  insight  into  the  errors  that  can  be  introduced  in 
the  determination  of  water  content  or  saturation  if  the  presence  of  layering  in  the  natural  system  is 
not  taken  into  account.  The  dielectric  constant-water  content  and  dielectric-water  saturation 
relationships  predicted  for  layered  systems  are  used  as  the  baselines  for  comparisons.  As  an 
example,  let  us  consider  a  system  (50%  fine  fraction)  with  thin  layering  perpendicular  to  the  wave 
propagation  direction.  In  Figure  3.2  are  shown  the  errors  that  would  result  if  water  content  and 
saturation  were  determined  from  dielectric  measurements,  and  the  system  was  assumed  to  be 
homogeneous;  i.e.  the  presence  of  layering  was  neglected.  The  errors  for  water  content  are  shown 
in  the  upper  plot,  while  the  errors  for  water  saturation  are  shown  in  the  lower  plot.  The  errors  for 
all  three  soil  systems  (sand-clay,  sand-silty  clay,  and  sand-silt  loam)  are  shown  for  use  of  both  the 
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Figure  3.1.  A  comparison  of  dielectric  constant  calculated  for  the  sand-clay  system  at  a  fine  fraction  of  0.50.  a) 
Dielectric  constant  versus  water  content,  b)  Dielectric  constant  versus  saturation. 


Difference  in  Water  Saturation 
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Figure  3.2.  a)  Difference  between  estimated  and  actual  water  content,  b)  Difference  between  estimated  and  actual 
water  saturation.  For  EMT  perpendicular  to  propagation  for  sand-fine  systems  at  a  0.50  fine  fraction. 


TP  model  and  the  Topp  equation.  As  can  be  seen  in  Figure  3.2,  the  errors  in  extracting  water 
content  can  be  as  high  as  0.064.  The  errors  in  extracting  water  saturation  are  even  greater  because 
of  the  effect  of  porosity  due  to  the  nature  of  the  binary  mixture  for  homogeneous  systems. 

Conclusions 

In  summary,  if  sand-fine  layers  exist,  using  either  the  TP  model  or  Topp  equation  to  extract  water 
content  or  saturation  can  produce  significant  errors.  Often  these  errors  fall  outside  the  acceptable 
limits  of  ±0.03  for  water  content  and  of  ±0.05  for  water  saturation.  For  all  three  layered 
scenarios,  the  TP  model  underestimates  water  content.  In  general,  the  Topp  equation 
underestimates  water  content  when  thin  layers  exist  parallel  to  the  EM  wave  propagation  but 
overestimates  water  content  when  either  thin  or  thick  layers  are  perpendicular  to  the  propagation 
direction.  However,  the  relationships  for  homogeneous  systems  tend  to  overestimate  water 
saturation  of  layered  systems.  The  errors  become  unacceptably  large  at  higher  water  saturations 
and  can  even  predict  water  saturations  that  have  impossible  values  which  are  greater  than  1. 

Numerical  Modeling  of  EM  Wave  Propagation  in  Layered  Systems 

Introduction 

In  Chan  and  Knight  (1998),  described  briefly  above,  we  model  the  relationship  between  dielectric 
constant  and  water  content/saturation  in  homogeneous  and  layered  systems  of  sand  and  fines.  In 

the  layered  models,  both  the  wavelength  (A)  of  the  EM  wave  and  the  average  thickness  (f)  of  the 
layers  are  accounted  for.  When  the  wavelength  is  much  larger  than  thickness  of  the  layers, 
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effective  medium  theory  (EMT)  is  used  to  calculate  the  average  dielectric  constant  of  the  layers. 
When  die  wavelength  is  much  smaller  than  the  thickness  of  the  layers,  ray  theory  is  used  to 
calculate  the  average.  The  region  over  which  each  of  these  relationships  is  valid  can  be  defined  in 
terms  of  the  ratio  k/t.  It  is  this  difference  between  the  EMT  and  ray  theory  relationship  that  is  the 
focus  of  this  study:  we  wish  to  determine  the  exact  region  (defined  by  X/t)  over  which  each 
relationship  is  valid.  While  it  is  clearly  important  to  account  for  the  presence  of  layering,  it  is  also 
necessary  to  determine  how  the  scale  of  the  layering  affects  wave  propagation,  and  the  resulting 
sampled  dielectric  constant 


Method  of  Numerical  Analysis 


The  variable  of  interest  in  this  study  is  X/t,  as  this  ratio  defines  the  transition  zone  between  a 
system  which  can  be  described  using  EMT  and  one  which  can  be  described  using  ray  theory.  The 
important  parameter  in  finding  the  transition  zone  between  EMT  and  ray  theory  is  X/t ;  i.e.  the  ratio 

between  the  wavelength  X  and  the  average  thickness  t  of  the  layers. 


Through  numerical  analysis  the  regimes  for  EMT  and  ray  theory  and  the  transition  zone  between 
them  for  EM  waves  were  determined  using  a  wave  propagation  matrix  method  written  by  Steven 
Cardimona.  This  program  although  based  on  wave  propagation  matrices  uses  the  Kennett  method 
for  stability.  The  program  calculates  the  travel  time  of  an  EM  wave  through  a  series  of  layers  with 
given  dielectric  constants  and  conductivities.  From  the  total  travel  time  and  the  total  thickness  of 
layers,  an  average  velocity  can  be  calculated  and  compared  to  the  velocities  derived  from  the 
dielectric  constants  predicted  from  EMT  and  ray  theory. 
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Wave  propagation  matrices  are  based  on  an  exact  solution  to  the  1-dimensional  wave  equation  for  a 
wave  traveling  through  a  series  of  layers.  As  described  by  Ward  and  Hohmann  (1994),  the  electric 
and  magnetic  fields  as  seen  in  Figure  3.3  can  be  represented  as  a  uniform  plane  wave: 

)+~Eie~ik' 5  (1) 

H  =  — (2) 

where  k(  is  the  complex  wave  number  in  the  i^1  layer, 

K  =  (3) 

co  is  the  angular  frequency,  HQ  is  the  permeability  of  free  space,  m  is  the  permeability  of  the  Ith 

layer,  £/  is  the  permittivity  of  the  Ith  layer,  z/  is  the  vertical  distance  to  the  bottom  of  the  j*  layer,  z 

is  any  vertical  distance  within  a  layer  at  which  the  field  is  measured,  +E  is  the  amplitude  of  the 
positive  traveling  electric  wave  in  the  i^1  layer,  and  ~E  is  the  amplitude  of  the  negative  traveling 
electric  wave  in  the  i^1  layer. 


Because  both  tangential  EM  fields  must  be  continuous  across  layer  boundaries,  the  EM  fields 
immediately  above  a  boundary  can  be  rewritten  in  terms  of  the  EM  fields  imnwtiafp  below  the 
boundary: 


E>«-1)=%  cosh  O’VO  -  Z.tf,  sinhO*,A) 

(4) 

cosh  (ikfo)  -  ^-Eyi  sinh(tt^) 

(5) 
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Figure  3.3  A  plane  EM  wave  normally  incident  to  a  series  of  layers.  Ey,  Hx,  and  it  are  a  right- 
handed  orthogonal  set. 


where  Zj  is  the  impedance  of  the  Ith  layer 


(6) 


and  hi  is  the  thickness  of  the  layer.  Since  equations  4  and  5  are  true  at  every  boundary  and 
since  they  can  be  written  in  matrix  form,  the  EM  fields  at  any  given  layer  is  simply  the  product  of 
all  the  layers  below  it: 


n 

cosh  (ik,hi) 

-Z#.sinh(/*A)‘ 

frr  1 

^>(4-1) 

/w 

=n 

i=i 

— —  sinh(ijfc^) 

.  3? 

cosh 

i - 

1  1 
d _ 

These  equations  have  been  programmed  to  calculate  the  average  dielectric  constant  of  an  EM  wave 
through  a  set  of  layers.  However,  one  must  proscribe  the  boundary  conditions  at  the  lowest 
boundary;  and  in  our  case,  we  have  set 


H,= 


m0 


(8) 


so  that  there  is  no  impedance  contrast  at  this  boundary  and  no  reflection  off  the  bottom  back  to  the 
surface  simulating  a  stack  of  layers  over  an  infinite  half-space  with  the  same  properties  as  the 
lowest  layer.  By  changing  both  the  thicknesses  of  the  layers  and/or  the  frequency  of  the  wave, 
wave  propagation  in  both  the  EMT  and  ray  theory  regimes  can  be  simulated. 


Data/Results 

Waveforms  were  numerically  propagated  at  50,  200,  and  750  MHz  through  layered  systems 
ranging  from  a  wavelength  to  layer-set  thickness  ratio  of  less  than  0.1  to  greater  than  10.  The 
dielectric  constants  of  the  alternating  equal  layers  were  15  and  5.  In  Figure  3.4  is  an  examples  of 
the  50  MHz  waveform  modeled  to  arrive  at  the  far  end  of  the  layered  system.  In  all  simulations  we 
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observe  travel  times  increase  as  the  layered  section  goes  from  a  few  thick  layers  (ray  theory  limit) 
to  many  thin  layers  (EMT  limit).  This  is  labeled  on  the  plot  in  terms  of  the  ratio  of  wavelength  1  to 
layer  thickness  t.  In  the  transition  zone  between  the  two  limits,  the  waveforms  become  very 
attenuated  due  to  scattering  and  the  travel  times  are  difficult  to  pick.  Scattering  also  causes  the 
wave  to  travel  more  slowly  than  theoretically  expected. 

In  Figure  3.5,  normalized  EM  velocity  versus  the  wavelength  to  layer  set  thickness  ratio  is  plotted 
for  50, 200,  and  750  MHz.  The  ray  theory  and  EMT  limits  are  shown  by  the  bold  lines.  Notice 
that  the  transition  zone  is  quite  narrow  and  is  slightly  higher  than  unity.  Other  calculations  indicate 
that  typical  soil  conductivities  affect  the  amplitude  of  the  waveforms  but  not  the  velocity  of  the 
wave.  Also,  the  transition  zone  moves  slightly  to  lower  values  when  the  proportion  of  high 
dielectric  material  is  more  than  half  and  moves  slightly  closer  to  one  when  the  proportion  of  high 
dielectric  material  is  less  than  half.  The  transition  zone  is  more  distinct  and  farther  from  one  when 
the  contrast  between  the  two  materials  is  small  but  is  more  diffuse  and  slightly  lower  when  the 
contrast  between  the  two  materials  is  large. 

Conclusions 

The  numerical  wave  propagation  modeling  confirms  that  the  EMT  relationship  or  the  ray  theory 
relationship  for  layered  materials  should  be  used  for  the  interpretation  of  field  data.  Because  the 
relationships  for  EMT  and  ray  theory  predict  different  dielectric  constants  for  the  proportions  of 
sand  and  fines,  one  must  determine  which  relationship  should  be  used.  The  EMT  relationship  is 
valid  when  the  wavelength  to  layer  set  thickness  ratio  is  large  while  the  ray  theory  relationship  is 
valid  when  the  ratio  is  small.  In  general,  the  transition  zone  between  the  two  limits  is  narrow  and 
occurs  just  below  unity. 
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Laboratory  Measurements 

Introduction/Experimental  Methods 

A  series  of  laboratory  experiments  were  conducted  with  systems  composed  of  alternating  layers  of 
sand  and  silt  These  experiments  were  designed  to  compare  the  theoretical  and  numerical 
predictions  (described  above)  to  results  obtained  through  measurements  on  geological  materials. 
We  were  interested  in  determining  whether  the  transition  from  effective  medium  theory  to  ray 
theory  could  be  observed  experimentally. 

The  measurement  apparatus  consists  of  a  center  rod  and  a  concentric  outer  shield.  The  center  rod 
is  1  cm  in  diameter  and  the  inner  radius  of  the  outer  shield  is  10  cm.  Both  the  center  rod  and  the 
outer  shield  are  segmented  into  10  cm  pieces.  The  center  rod  segments  screw  together  while  the 
outer  shield  segments  stack  on  top  of  each  other.  The  outer  shield  segments  are  also  scored  every 
half  centimeter  and  are  sealed  with  O-rings.  The  total  height  of  the  apparatus  is  1.20  m.  However, 
because  the  apparatus  is  segmented,  the  total  height  does  not  need  to  be  used  for  every 
measurement.  The  bottom  segment  is  fixed  to  a  base.  In  the  base  is  a  coaxial  cable  connection 
which  connects  the  coaxial  cell  to  the  data  collection  instrument  Also  in  the  base  is  a  mechanical 
switch  which  causes  a  short  between  the  center  rod  and  the  outer  shield.  The  base  also  has  a  valve 
to  allow  fluid  flow. 

To  pack  the  cylinder,  a  premeasured  weight  of  soil  is  placed  in  the  bottom  segment  of  the  cylinder 

and  tamped  down  until  it  fills  the  desired  layer  thickness.  For  saturated  layers,  the  calculated 

amount  of  water  is  poured  onto  a  mesh  which  lies  on  top  of  the  soil  and  distributes  the  water 

evenly  throughout  the  soil.  This  process  is  repeated  until  the  desired  number  of  layers  is  reached. 

* 

Cylinder  segments  are  added  to  the  apparatus  as  they  are  needed. 
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Data  Collection  and  Results 

The  data  were  collected  using  a  Tektronix  1502C  metallic  cable  tester  TDR.  This  instrument  is 
controlled  by  a  486  PC  using  an  RS232  interface  and  the  program  TDR-Main  written  by  David 
Redman  at  the  University  of  Waterloo.  This  program  records  the  amplitude  of  the  TDR  trace 
between  two  specified  times.  The  data  arc  converted  to  a  text  format  and  imported  to  MatLab  for 
data  analysis. 

For  each  measurement  three  TDR  traces  are  collected — an  open  trace,  a  trace  where  the  cable  is 
shorted  at  the  base  of  the  cylinder  using  the  mechanical  switch  mentioned  above,  and  a  trace  where 
the  center  rod  and  the  outer  sheath  are  shorted  using  a  metal  rod  immediately  above  the  top  layer  of 
soil.  The  three  traces  are  plotted  together.  The  times  when  the  shorted  traces  deviate  from  the 
open  trace  are  used  to  calculate  the  travel  time  through  the  sample  and  in  turn  the  average  dielectric 
constant. 

In  order  to  determine  the  accuracy  of  our  apparatus,  measurements  were  taken  on  air  and  water. 
Measurements  were  taken  at  every  layer  (ie  every  10  cm)  of  the  cylinder.  The  results  are  within  a 
few  percent  of  the  expected  values  of  approximately  1  and  80,  respectively. 


Data  were  collected  on  alternating  layers  of  dry  sand  and  saturated  silt.  Measurements  were  taken 
on  layers  20  cm,  10  cm,  and  1  cm  in  thickness.  Assuming  a  frequency  of  750  MHz,  this 
corresponds  to  a  ratio  of  wavelength  to  layer-set  thickness  (A/t)  of  approximately  0.3,  0.6,  and 
6.0,  respectively.  The  dielectric  constant  of  each  layer  was  calculated  using  the  TP  model.  The 
resulting  values  were  then  used  to  calculate  the  theoretical  EMT  and  ray  theory  predictions.  Figure 
3.6  shows  normalized  EM  wave  velocity  versus  A/t.  This  plot  is  very  similar  to  Figure  3.5  and 
confirms  that  both  EMT  and  ray  theory  limits  are  measurable.  However,  more  data  are  needed  to 
fill  in  the  transition  zone  more  fully. 


(^«H-XW3)/(P*JnsBaW-XIVa) 
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Figure  3.6.  Normalized  velocity  versus  wavelength  to  layer  set  thickness  ratio  from  experimental  study.  Half  the  layers  are  wet  silt  and  half 
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Conclusions 

Laboratory  data  confirm  that  the  EMT  relationship  or  the  ray  theory  relationship  for  layered 
materials  should  be  used  for  the  interpretation  of  dielectric  constant  We  observe  a  transition  from 
EM  velocities  predicted  by  ray  theory,  at  low  values  of  X/t,  to  EM  velocities  predicted  by  effective 
medium  theory,  at  large  values  of  X/t.  The  transition  zone  between  these  two  limits  was  found 
close  to  1,  as  predicted  by  numerical  modeling. 

CONCLUSIONS 

Measurement  of  the  dielectric  constant  of  geological  materials  can  be  used  to  estimate  hydraulic 
properties  -  but  the  accuracy  of  the  determined  hydraulic  properties  depends  on  the  use  of  valid 
relationships  to  link  dielectric  properties  to  hydraulic.  Most  studies  to  date  have  been  based  on 
laboratory  and  theoretical  results  that  use  or  assume  homogeneous  samples.  In  real  materials  there 
are  likely  to  be  layers  of  different  materials  present;  these  will  affect  both  the  measured  dielectric 
constant  and  the  relationship  to  hydraulic  properties. 

Our  theoretical  studies  showed  that  significant  error  can  be  incurred  if  relationships  for 
homogeneous  systems  such  as  the  TP  model  or  Topp  equation  are  used  to  extract  water  content  or 
saturation  from  the  dielectric  constant  measurements  of  layered  systems.  Additionally,  numerical 
EM  wave  propagation  modeling  showed  that  the  wavelength  to  layer  set  thickness  ratio  of  layered 
systems  can  also  influence  the  average  dielectric  constant  of  a  layered  system.  Finally,  laboratory 
studies  confirmed  that  both  the  EMT  limit  and  the  ray  theory  limit  can  be  measured  in  sand  and 

silt  The  numerical  and  experimental  results  show  a  transition  from  ray  theory  to  EMT  at  X/t  close 

to  1.  We  conclude  that  heterogeneity  such  as  layering  must  be  accounted  for  in  order  to  accurately 
interpret  dielectric  measurements. 
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4 

Electrical  Conductivity  Estimates  from  GPR  Data 

A  central  theme  in  our  research  is  the  desire  to  use  dielectric  information,  obtained  from  GPR  data, 
to  determine  hydrogeologic  properties.  It  is  clear  however  that  there  is  not  a  simple  transform 
between  these  two  types  of  properties.  Additional  information  is  required.  The  idea  that  we  have 
investigated  is  that  there  exists  in  the  GPR  data  itself  additional  information  that  is  presendy  under¬ 
utilized.  Specifically,  GPR  data  contain  information  not  only  about  dielectric  properties,  but  also 
about  electrical  conductivity.  This  work  is  described  in  detail  in  the  paper  by  Rea  and  Knight 
(1997),  submitted  to  Geophysics. 

The  observed  attenuation  in  GPR  data  contains  information  about  the  electrical  conductivity  of  the 
subsurface.  By  quantifying  all  contributions  to  the  attenuation  recorded  in  the  GPR  data  we  can 
isolate  and  estimate  the  electrical  conductivity.  By  far  the  most  ill-constrained  parameter  in  this 
method  is  the  estimation  of  power  loss  associated  with  the  partitioning  of  energy  at  interfaces  -  as 
some  energy  is  transmitted  and  some  is  reflected  at  each  interface  across  which  there  is  a  change  in 
dielectric  properties.  We  have  concluded  that  the  best  approach  is  to  attribute  all  loss  due  to  energy 
splitting  at  interfaces  to  conductive  loss;  our  analysis  of  the  attenuation  will  then  result  in  our 
obtaining  a  maximum  or  upper  bound  for  the  electrical  conductivity. 

We  have  found  that  this  method  can  be  used  to  quantify  both  lateral  and  vertical  changes  in 
electrical  conductivity.  We  first  tested  the  method  at  a  site  where  there  exists  lateral  variation  in 
electrical  conductivity  (a  sand  and  gravel  unit  adjacent  to  a  day-rich  unit).  At  this  location  we  had 
well  control  and  conducted  a  DC  resistivity  survey  along  the  same  line.  Conductivity 
measurements  were  also  available  from  Time  Domain  Electromagnetic  soundings  at  locations 
within  100  m  of  the  line.  The  results  that  were  obtained  from  our  analysis  of  the  100  and  200 
MHz  GPR  data  are  listed  in  Table  1  and  compared  with  those  from  the  DC  resistivity  and  the 
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Table  1  :  Conductivity  results  from  GPR  decay  rate  analysis,  and  DC  resistivity  and 
TDEM  surveys. 


sedimentary 

unit 

100  MHz  GPR 
o  (mS/m) 

200  MHz  GPR 
a  (mS/m) 

DC  resistivity 
a  (mS/m) 

TDEM 
a  (mS/m) 

3.5 

4.3 

0.2  -  1.6 

0.1  -  5.0 

27 

51 

9.5  -  47 

25-50 

TDEM  surveys.  As  can  be  seen,  we  obtained  electrical  conductivity  estimates  from  the  GPR  data 
that  are  in  very  good  agreement  with  those  obtained  from  other  methods. 

At  the  second  test  site,  a  sand  and  gravel  unit  overlies  a  massive  clay,  producing  a  large  vertical 
variation  in  electrical  conductivity.  We  were  able  to  see  the  response  of  each  of  these  units  and 
determined  conductivity  values  that  agreed  well  with  other  measurements  at  the  site. 

The  attenuation  seen  in  GPR  data  has  long  been  interpreted  in  a  qualitative  way  by  practitioners  in 
terms  of  conductivity;  i.e.  the  observation  of  reduced  signal  penetration  indicates  a  more 
conductive  region.  What  we  have  shown  in  this  study  is  that  it  is  possible  to  extract  a  reasonable 
estimate  of  electrical  conductivity  from  an  analysis  of  the  observed  attenuation.  This  provides  us 
with  valuable  complementary  information  and  an  additional  parameter  to  help  constrain  the  problem 
of  obtaining  a  reliable  model  of  the  material  or  hydrogeologic  properties  of  the  subsurface  from 
GPR  data. 
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5 

Conclusions  and  Recommendations 

The  goal  of  our  research  is  to  develop  an  understanding  of  the  link  between  dielectric  properties 
and  hydraulic  properties  of  the  subsurface.  Measurements  of  dielectric  properties  of  the  subsurface 
can  be  obtained  non-invasively  from  ground  penetrating  radar  data.  From  an  applied  perspective 
our  research  can  be  described  as  addressing  the  question:  How  can  we  use  these  measurements  of 
dielectric  constant  to  characterize  hydraulic  conductivity  at  a  site? 

Our  approach  in  this  three-year  research  project  has  been  a  combination  of  field  experiments, 
laboratory  experiments,  and  theoretical  modeling.  Regardless  of  the  approach,  the  emphasis  has 
been  to  work  with  well-characterized  materials,  or  field  sites,  so  that  we  can  study  in  detail  the 
relationship  between  the  hydrogeologic  properties  of  geological  materials  and  their  dielectric 
properties. 

There  are  two  types  of  hydrogeologic  information  that  have  been  the  focus  of  our  research.  One 
type  of  information  concerns  obtaining  estimates  of  the  actual  magnitude  of  hydrogeologic 
parameters;  e.g.  what  is  the  porosity  of  this  zone?  what  is  the  water  content  of  the  top  meter?  ITie 
other  type  of  information  concerns  the  spatial  variability  in  hydrogeologic  properties;  e.g.  over 
what  distance  do  we  expect  to  see  continuity  in  permeable  zones  -  meters  or  10’s  of  meters? 

The  first  type  of  information  can  be  obtained  by  taking  measurements  of  dielectric  properties  of  a 
volume  of  the  subsurface  and  transforming  it  to  obtain  estimates  of  the  hydraulic  properties  of  that 
volume.  The  accuracy  with  which  this  can  be  done  obviously  depends  heavily  on  the  validity  of 
the  relationship  that  is  used  to  transform  dielectric  measurements  to  hydraulic  properties.  Very 
commonly  the  assumption  is  made  that  the  sampled  volume  is  homogeneous;  this  being  the 
underlying  assumption  is  most  relationships  used  to  link  dielectric  and  hydraulic  properties.  Our 
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theoretical,  numerical  and  laboratory  studies  have  addressed  the  effect  of  the  presence  of  layering, 
and  the  importance  of  considering  both  the  thickness  of  the  layers  t  and  the  wavelength  of  the 
measurement  X.  We  find  that  when  X/t  is  less  than  1  ray  theory  should  be  used  is  the  interpretation 
of  dielectric  properties  for  hydraulic  properties  and  when  X  /t  is  great  than  1,  an  effective  medium 
theory  should  be  used.  We  have  calculated  the  error  in  involved  in  neglecting  the  effect  of  layering 
and  find  that  it  can  be  significant  The  results  from  this  aspect  of  our  research  suggest  that 
heterogeneity  should  be  considered  in  interpreting  dielectric  data;  even  if  this  only  involves 
obtaining  estimates  of  the  likely  errors  involved  in  the  interpretation  of  the  data  for  hydrogeologic 
properties.  In  the  natural  world  the  assumption  of  heterogeneity  is  much  more  likely  to  be  valid 
that  the  assumption  of  homogeneity. 

The  second  type  of  required  information  involves  describing  the  spatial  heterogeneity  of  the 
subsurface.  Often  at  a  site  there  are  some  direct  measurements  of  subsurface  properties,  but 
“filling  in”  between  these  measured  data  points  requires  knowledge  of  the  structure  or  spatial 
variability  of  the  subsurface,  at  that  specific  site.  We  have  found  that  geostatistical  analysis  of 
GPR  data  can  be  used  to  characterize  spatial  variability  in  the  subsurface,  and  have  had  tremendous 
success  in  analyzing  data  from  a  variety  of  depositional  environments.  This  is  an  approach  that  we 
highly  recommend  being  used  at  sites  either  in  the  initial  planning  stages,  or  to  assist  in  the 
integration  and  interpretation  of  data.  If  it  is  possible  to  collect  a  GPR  image  at  a  site,  that  image  is 
an  excellent  source  of  information  about  the  heterogeneity  of  the  subsurface.  Further  research  in 
this  area  is  needed  to  determine  more  specifically  the  link  between  the  GPR  image  and  the 
hydrogeologic  properties  -  i.e.  are  we  primarily  imaging  variation  in  water  content,  clay  content, 
permeability,  etc.?  This  is  the  critical  question  that  must  be  answered  in  order  to  further  the 
effective  use  of  GPR  for  obtaining  information  about  the  nature  and  distribution  of  subsurface 
materials. 
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The  combination  of  our  field,  laboratory  and  theoretical  work  has  provided  us  with  new  ideas,  and 
new  approaches,  to  obtaining  hydrogeologic  information  from  dielectric  measurements.  In  many 
applications  in  the  earth  sciences  and  engineering  we  are  faced  with  the  challenge  of  characterizing 
a  complex,  highly  heterogeneous  system.  As  we  further  our  understanding  of  the  link  between 
dielectric  properties  and  hydrogeological  properties,  we  improve  our  ability  to  use  dielectric 
measurements  as  a  non-invasive  means  of  characterizing  the  scale-dependent  hydraulic  properties 
of  the  subsurface. 
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